AUTOMATED SYSTEM FOR ANALYZING THE ACTIVITY OF INDIVIDUAL NEURONS 
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ABSTRACT 


This paper presents a signal processing system that i) provides an efficient and reliable instrument for 
investigating the activity of neuronal assemblies in the brain and ii) demonstrates the feasibility of generating the 
rnmmaiul signals of prostheses using the activity of relevant neurons in disabled subjects. The system operates on- 
line in a fully automated manner and can recognize the transient waveforms of several neurons in extracellular 
neurophysiological recordings. Optimal algorithms for detection, classification, and resolution of overlapping 
waveforms are developed and evaluated. Full automation is made possible by an algorithm that can set appropriate 
decision thresholds and an algorithm that can generate templates on-line. The system is implemented with a fast 
IBM PC compatible processor board that allows on-line operation. 


INTRODUCTION 


Two equally important reasons have recently increased considerably the significance of processing the signals 
generated by neurons. 

1) The expanding applicability of neural networks to diverse engineering problems has raised the interest in 
neurophysiological investigations that aim to study the collective behavior of neuronal assemblies. It seems clear 
that advances in areas such as pattern recognition, memory storage, speech processing, computer vision, and 
control will be possible by a better understanding of biological neural systems, specially the human brain. Trus 
promise has been recognized in the Congressional resolution designating the 1990 s the Decade of the Brain. 
Furthermore the National Academy of Sciences indicated that neuroscience stands at the threshold of a significant 
expansion. However, as neuroscientists in The Johns Hopkins Medical School and other leading universities agree, 
an important prerequisite is instrumentation for observing the electrical activity of individual neurons and their 
assemblies. 

2) In aiding the disabled, various kinds of prostheses have been effective and have led to increased attention to the 
field. In current prostheses, the motion command for a joint assisted or replaced by a prosthesis is generated by 
another joint: for example, the command for the legs of a paraplegic is initiated with switches operated manual y 
by the subject. A recent direction of research, named neural prostheses, aims to generate the commands in a more 
convenient, natural and potentially more effective manner. In a neural prosthesis, the goal is to obtain the 
command directly from neurons or muscle fibres that are relevant to the target joint. If this can be achieved by 
obtaining the commands from the very neurons that once controlled the disabled joint, the most direct and natural 
link between the motion intent and the target joint can be developed. Here again, the prerequisite is 
instrumentation for obtaining the activity of neurons. 


Sensory or motor information is processed by biological systems in the form of a distributedrepresentation 
supported by a large number of neurons interconnected with excitatory and inhibitory synapses. The functional 
activity of an individual neuron depends on the strength of the synapses that provide excitatory or inhibitory input 
from other neurons, and on the activity of these neurons. When the net input is excitatory , above a threshold level, 
the neuron fires and generates an action potential across its membrane, that lasts about 1 ms. Ongoing activity in a 
neuron is manifested by a sequence of action potentials that can be recorded with an extracellular electrode. The 
main advantage over an intracellular electrode is the ability to record from more than one neuron at the same 
time but the extracellular electrode also allows recording without damaging the neurons. The cost of these 
benefits is the requirement for sorting the interleaved neural spike trains (Fig. 1) to determine the firing instants 
of individual neurons. 
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Fig. 1. Recordings at two different noise levels. The number on each spike indicates the 
corresponding neuron. The SNR of a spike in this study is defined as the rms value of the spike in 
units of standard deviation of noise. The SNR of the smallest spike (2) is about 4 in the top trace 
and about 2 in the bottom trace. Superpositions of pairs of spikes are labeled as 3-1 on the top 
trace and 2-1 on the bottom trace. The sampling rate was 32 KHz. 


Different neurons generate distinct spike waveforms in the recording, due to differences in their dendntic 
geometry and the impedances of the medium connecting them to the electrode. The activity of individual neurons 
can be determined by sorting the different types of neural waveforms. An essential challenge in extracellular 
recordings is the relatively low si gnal-to- noise ratios (SNR) that can occur in many cases. The background noise is 
mainl y due to the activity of a large number of distant neurons resulting in a considerable overlap between the 
spectra of waveforms of interest and noise. Furthermore, since the noise process is primarily made of the 
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accumulation of low amplitude spikes, the autocorrelation function of noise has significantly high coefficients at 
lags as large as the average duration of a neural spike (about 1 ms). Moreover, two neurons often fire 
concurrently, their spikes overlap in the recording (Fig. 1) , causing additional difficulty in recognizing the 
activity of individual neurons. The neural waveform recognition problem is a typical example of detection and 
classification of transient patterns embedded in colored noise. 

In recordings for brain research, it is preferable to achieve neural spike sorting on-line because immediate 
feedback on the recording allows better control of experimental conditions and recording quality, reducing the 
time requirement for the neuroscientist and the animal subject In neural prosthesis applications on-line operation 
is essential. This requires a signal processing system that can operate in a fully automated manner, without human 
supervision. We developed signal processing and pattern recognition algorithms as well as a hardware 
implementation that operate with theoretically optimal recognition performance, on-line, in a fully automated 
manner [1-5]. 

The data used in the development of this analyzer were recorded from the primate cortex using a filter pass- 
band of 10 Hz to 10 KHz, a 12 bit A/D converter and a sampling rate of 32 KHz. 


ALGORITHMS 

In order to determine the activity of individual units, an automated system should perform three main tasks: i) 
discrimination of waveforms from noise (detection), ii) discrimination between waveforms of different units 
(classification), and iii) separation of overlapping waveforms (resolution of superpositions). 


Detection 


The analyzer performs detection by computing the power p(k) of the signal N(k) within a running window of 
length/I : 

n-1 

p(k) = ZN(k-i) . (1) 

i=0 


A waveform is detected when the power exceeds an appropriate acceptance threshold which is set with an 
automated algorithm. Power detection yielded considerably better results than the widely used amplitude detection 
in tests comparing the detection performance as a function of SNR. The SNR was defined as the ratio of the signal 
ims value (computed with n samples) to the standard deviation of noise in the record. With the detection 
thresholds set for a false positive rate of less than 0.1%, power detection was 100% and 94% correct at SNRs of 3 
and 2, while amplitude detection provided 95% and 71% correct detection respectively. The system generates, on- 
line, a template for each unit to be used in classification. When the templates become available, the detection of the 
corresponding waveform types can be improved with matched filtering. The improvement obtained with matched 
filtering can be as high as 40% more correct detection than the power detection technique, at low SNR levels. 


flaggifinttirw 

Several methods for neural spike classification, ranging from amplitude discrimination to principal 
components and minimum mean-square-error, have been suggested [6-13]. In previous comparison studies [7] 
principal components and template matching using Euclidean distance were found to be the best for spike sorting 
in noisy data. Because the statistically optimal Bayesian classification can be achieved by template matc hin g and 
because current technology allows its on-line implementation, we focused on neural spike sorting by template 
matching. In view of the diversity of applicable classification techniques, the emphasis of this project was placed 
on developing the theoretically optimal approach that could provide minimal noise sensitivity. Therefore, the 
optimal template matching approach was investigated and evaluated. 
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In the template matching method, each waveform is represented by n consecutive samples digitized on the 
waveform and it can be viewed as an n -dimensional vector. The waveform vector of a given neuron will vary 
somewhat each time that this neuron fires, due to additive random noise. The waveforms of the same neuron will 
form a cluster in sample space, around the centroid that would be observed in the absence of noise. When the 
Hiotrihi.tifm of noise amplitude is Gaussian, as in neural recordings, the optimal Bayesian classification can be 
achieved by computing the mean of the cluster of each neuron and by setting a decision boundary around each 
mean with a distance metric that depends on the covariance matrix of noise. 

The probability density p(x) of a multivariate Gaussian cluster distribution is. 


1 

p(x) expKx-m^C'Vx-myZ], 

( 2«) n/2 ICI iy2 


( 2 ) 


where x is the n -dimensional random waveform vector, C is the covariance matrix of x vectors in that cluster, ICI 
is the determinant of C, and m is the mean vector in the cluster. In the multiclass problem, each class has its own 
cluster and probability density. 


Let the data have K different classes represented by Wj with the class index i = 1.....K. Multiclass Bayesian 
classification is performed with discriminant functions that are based on the class densities and a priori 
probabilities of the classes. A convenient choice of discriminant function is: 


g.(x) = logljKxlw.) + logP(w.), 


(3) 


where, g (x) is the discriminant function, p(xlwj) is the probability density of class i, and P(wj) is the a priori 

orobability of class i. The class to which a candidate pattern belongs is determined by computing the values of the 
discri minant function using the pattern’s vector x for each class. The pattern is assigned to the class with the 
highest discriminant function value. 


When each class has a multivariate Gaussian distribution, the expression fix gj(x) becomes. 
g.(x) = -(x-m. ) T C _1 . (x-m .)/2 - (nlog2jc)/2 - log)C.I)/2 + logPCw.). 


(4) 


This expression can be further simplified because the (nlog2jt)/2 term is the same for each class, the covamnce 
matrix C is the same for all classes and the a priori probability of each class is assumed to be the same. Therefore, 

classification with this discriminant function becomes equivalent to classification according to the minimal value of 
the quantity: 

d i 2 = (x-m i )‘C* 1 (x-m i ) < 5 ) 

which is known as the squared Mahalanobis distance between x and nr. Constant Mahalanobis distance contours 

are ellipses centered around the mean of each class; these ellipses coincide to equal density contours on the 
multivariate distributions. In most applications, it is desirable to have the option of leaving some patterns 
unclassified or rejecting them. To do so. only patterns that have a distance below an acceptance threshold are 
classified. This is equivalent to setting an elliptical decision boundary around the mean of each class. 

None of the reported neural spike sorting applications have used the Mahalanobis distance in classification 
because of its computational burden. The Euclidean distance which has been commonly used, is equivalent to 
setting a circular decision boundary and provides suboptimal results. The extent of performance loss due to the use 
of F uc |i d'*” n distance depends on the covariance matrix of noise, the noise level of the data and the similarity 
between different classes.The Euclidean distance can provide satisfactory results regardless <rf the covariance 
matrix of noise if the noise level is relatively low and the clusters of different classes are sufficiently apart. Our 
ISSiS 32-ampte .empires and 5 different neural apike claaaea embedded in typreal neunl re emtog 
noise, showed that if the Euclidean distance between the means of the two closest clusters is more than 14 standard 
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deviations of noise, perfect classification can be obtained But as the clusters get closer and/or the noise level 
increases, the performance drops and the loss, referred to the optimal case, can reach up to 25%. 

The extent to which Euclidean distance deteriorates the classification performance depends on the covariance 
matrix. The higher the autocorrelation in noise, the higher the eccentricity of the elliptical distributions and the 
lower the performance will be with Euclidean distance. However, if the noise had no autocorrelation, the 
covariance matrix would be diagonal and the Euclidean distance would provide optimal c la ssification. Therefore, 
the system that we developed whitens the data before classifying the waveforms. Whitening is achieved with a 
digital FIR or IIR filter whose coefficients are determined by modeling the noise in the recording. The model is 
an autoregressive moving average (ARMA) model [13] and its inverse provides the whitening filter. 

The contribution of whitening to template matching performance was evaluated using a test data set with 5 
different spike types and prior knowledge of exact templates. With the acceptance threshold set to provide fewer 
than 0.1% false positives, the correct classification performance on raw data was 96% and 74% correct, at SNRs 
of 2 and 1 respectively, while on whitened data the performance increased to 100% and 91% correct respectively. 


Resolution of superpositions 

When two neurons fire simultaneously their waveforms overlap and generate a complex waveform that does 
not match any of die templates of individual neuron waveforms. Such superpositions occur with considerable 
frequency depending on the number of neurons in the recording, their firing rates, the duration of the action 
potentials, and the timing relations between the action potentials of individual neurons. Failure to recognize the 
spikes that overlap can cause underestimation of the firing rates, and can affect severely the analytic measures of 
interevent timing, such as cross-correlation between neurons or inter-spike interval histograms. In a typical 
recording containing waveforms from 5 neurons firing at moderately high rates, about one third of the 
waveforms may overlap. 

The neural waveform analyzer that we developed includes a superposition resolution algorithm that essentially 
subtracts each template from the complex waveform and attempts to classify the remainder. The resolution is 
performed on whitened data and provides an optimal template matching approach for this task. The performance 
of this algorithm was evaluated in tests where all templates were available. With the acceptance threshold set to 
provide fewer than 0.1% false positives, the correct resolution performance was 100% and 95% correct at SNRs 
of 3 and 1.5 respectively, on whitened data. 


Besides the three mam functions of detection, classification and superposition resolution, two other functions 
are needed for a fully automated on-line system. The first one is automated threshold setting that requires a 
segment of only noise in the record. The decision thresholds for the power value in detection and the Euclidean 
distance in classification have to be set in accordance to the level of noise in the recording. Therefore, in order to 
set the thresholds automatically, a section that contains wily noise has to be segmented from the recording. 
Furthermore, this has to be done at the beginning of the process, without using any of the detection nv-hniqi ^» s 
mentioned above because thresholds are not yet available. We developed an iterative algorithm that can separate 
the noise from all other transient waveforms, without using templates or decision thresholds. This algorithm is 
based on the fact that the amplitude of the noise process is normally distributed. The noise segmentation algorithm 
provided a dequate results in test data sets that had SNR levels ranging from 1 to 10 and total spike rates of 5 to 
160 spikes per second. 


Template generation 

The second function required for full automation is template generation. The system should be able to observe 
the incoming data and determine a template for the waveform of each neuron. In the system that we developed, 
this is achieved with unsupervised clustering techniques that provide two different approaches: sequential or 
simultaneous. In the sequential clustering approach, the first spike detected becomes the first template, 
representing the waveforms of the first neuron (type 1). The second spike is compared to this template; if the 
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distance is lower than the acceptance threshold, the second spike is classified as type 1 and the template is updated 
bv averaging If the distance is greater than the acceptance threshold, the second spike waveform ts used as a 
KSplatTrepresenting^teSnd type of spike waveform. Each subsequent spike is 
either die corresponding template is updated or a new type is initiated. The sequential clustering algorithm 
provided appropriate templates in the data sets that we used for evaluating the system. We are currently 
investigating simultaneous clustering algorithms that have potential for increasing theclustenng performance trt 
very low SNR levels. Simultaneous clustering algorithms use clustering criteria applied simultaneously to a larg 
number of waveforms obtained at the start of the recording. The templates generated by such algorithms do not 
depend on the detection sequence of waveforms and the effects of noise at low SNR levels are reduced. 


Complete system 

The block diagram of the complete system is shown in Fig. 2. At the beginning of the recording, an initial 
section of the digitized data is processed by the noise segmentation algorithm that provides segments of only noise 
to the system. The thresholds for detection, classification and superposition resolution are set automatically using 
the extracted noise segments. The noise segments are also used by the ARMA modeling algorithm which intum 
provides the coefficients of the whitening filter. After these preprocessing steps that last several seconds, the 
remaining recognition functions are performed on-line using a double buffering input arrangement. 



Fig. 2. Block diagram of complete system. Shaded components represent preprocessing 
algorithms. 


At the start of the recognition process, since templates are not available, waveforms are detected with the 
power detector that passes indicates the occurrence times of the detected events to the classifier. The classifier 
applies sequential clustering to the whitened data to generate and update templates for each neuron in the 
recording and uses the templates for classifying the waveforms on the whitened data. Each time that a wavrform 
matches a template, the corresponding type and occurrence time is recorded and the matching template is updated. 
If the waveform does not match any template, it is first assumed to be a superposition and placed in the 
corresponding stack. 


The superposition resolution algorithm attempts to resolve each waveform in this stack by using the templates. 
If the unknown waveform can be resolved in terms of two templates, these two templates are updated and two 
spikes are recorded. If the waveform cannot be resolved then it is considered to represent a new neuron that 
started to fire and this waveform becomes a new template. 
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In some cases, waveforms that cannot be classified or resolved can be artifacts or spikes of very inactive 
neurons that are not worth pursuing. The template controller monitors the activity of each type by computing the 
number of times each template has matched a waveform in the recent past. If the activity of a given template is 
lower than a preset level, that template is eliminated. 

When the template for a neuron is available, the detection of that type is performed by matched filtering, 
while power detection is kept active in order to detect new types as well as superpositions. 


HARDWARE 

The system is implemented on an IBM PC compatible, floating-point processor board developed in The Johns 
Hopkins Applied Physics Laboratory. This implementation allows about 40 MFLOPS operation for most of the 
functions of the system such power detection, whitening, classification, superposition resolution, matched filtering 
and template update, using assembly language. Each of the two 4 MBytes memory banks of the processor is 
connected to a parallel I/O port that can transfer data at a rate of 80 MBytes per second. The recorded signal is 
digitized with a commercial A/D board that stores the data directly on one memory bank of the processor using 
DMA through one of the parallel I/O ports, without passing from the IBM PC and without requiring time from 
the CPU of the processor. The on-line recognition results obtained by the processor are passed to the IBM PC 
through the second parallel I/O port for display and archiving. The human interface provides displays of the raw 
data, whitened data, running power values, matched filter outputs, templates, clusters projected on two 
dimensions, as well as measures of the data quality and the recognition difficulty. Further display functions such as 
correlation histograms, interval histograms, and raster plots of spikes can also be implemented. 


CONCLUSION 

We developed a fully automated system that can recognize the transient waveforms generated by several 
neurons in an extracellular recording. The most significant contributions of this system are i) theoretically optimal 
operation that provides minimal noise sensitivity, ii) the ability to generate all operational parameters (e.g. 
templates and thresholds) automatically and on-line, allowing its use with minimal human supervision, and iii) 
resolution of superpositions, providing an indispensable tool for more complete data acquisition. 

Furthermore, since the processor is a general purpose computation tool, its use is not limited to only sorting 
the waveforms. After the recording, the same hardware programmed with a high level language such as C, can be 
used for investigating the collective behavior of many neurons. By allowing both waveform recognition and 
further neurophysiological investigation, this system provides a cost effective instrument for neuroscience 
research. 

The fully automated and on-line operation is a unique property of this system that shows the feasibility of 
reliable on-line recognition of neural activity for neural prosthesis applications. In neuroscience research 
applications, fully automated on-line operation enables more efficient recording and processing. This is the result 
of the reduced human supervision and the immediate feedback that the system can provide to the user. By 
eliminating the need for the constant supervision that available systems require, the system that we developed 
allows more time and focus to the neuroscientist for the proper management of the experiment Moreover, by 
recognizing and reporting immediately the activity of neurons, the system provides valuable feedback and 
guidance for the recording. Research decisions that are made several days after the recording, can be made during 
the recording owing to the immediate, automated results. This more efficient operation can lead to a reduction in 
the use of laboratory animals. 

This system is an example of the transfer of military technology to civilian and commercial applications. 
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